You can complete all the coding exercises in the JAGS Primer using simple R scripts. However, you might be wondering about how we created the JAGS Primer key you are looking at right now. We did this using Yihui Xie’sknitr package in R. This can be a highly useful tools for organizing your Bayesian analyses. Within the same document, you can:
Describe the model and specify the joint distribution using \(\LaTeX\) and RMarkdown.
Write the R and JAGS code for implementing your model in JAGS.
Run the model directly in JAGS.
Summarize the convergence diagnostics and present results.
Add anything else pertinent to your analysis.
Best of all, knitr can produce beautiful html files as output, which can be easily shared with collaborators. We encourage you to become familiar with knitr. We recommend Karl Broman’s knitr in a nutshell as an excellent introductory tutorial.
Exercise 1: Factoring
There is no x in the posterior distribution in equation 4. What are assuming if x is absent? Draw the Bayesian network, or DAG, for this model. Use the chain rule to fully factor the joint distribution into sensible parts then simplify by assuming that \(r\), \(K\) and \(\tau\) are independent.
There is no \(x\) because we are assuming it is measured without error.
Fig 1. Bayesian network for logistic model.
Factoring the joint using the chain rule: \[[r,K,\tau, \mathbf{y}]=[y\mid r,K,\tau][r|K,\tau][K|\tau][\tau]\]
Simplifying under the assumption of \(r,K,\tau\) are independent: \[[r,K,\tau, \mathbf{y}]=[y\mid r,K,\tau][r][K][\tau]\]
We would obtain the same result by inspecting the DAG and using the rules:
All random variables on at the head of arrows must be on the lhs of a conditioning symbol.
All random variables at the tails of arrows must be on the rhs of a conditioning symbol.
All random variables at the tails of arrows with no arrows coming into them must be expressed unconditionally.
Exercise 2: Can you improve these priors?
A recurring theme in this course will be to use priors that are informative whenever possible. The gamma priors in equation 4 of the JAGS Primer include the entire number line \(>0\). Don’t we know more about population biology than that? Lets, say for now that we are modeling the population dynamics of a large mammal. How might you go about making the priors on population parameters more informative?
A great source for priors in biology are allometric scaling relationships that predict all kinds of biological quantities based on the mass organisms (Peters, 1983; Pennycuick,1992). If you know the approximate mass of the animal, you can compute broad but nonetheless informative priors on \(r\) and \(K\). This might leave out the social scientists, but I would trust the scaling of \(r\) for people if not \(K\).
In the absence of some sort of scholarly way to find priors, we can at least constrain them somewhat. There is no way that a large mammal can have an intrinsic rate of increase exceeding 1 – many values for \(r\) within gamma(.001, .001) are far large than than that and hence are complete nonsense. We know \(r\) must be positive and we can put a plausible bound on its upper limit. The only requirement for a vague prior is that its “\(\ldots\) range of uncertainty should be clearly wider that the range of reasonable values of the parameter\(\ldots\)” (Gelman and Hill, 2009, page 355), so we could use \(r\) ~ uniform(0, 2) and be sure that it would be minimally informative. Similarly, we could use experience and knowledge to put some reasonable bounds on \(K\) and even \(\sigma\), which we can use to calculate \(\tau\) as \(\tau=\frac{1}{\sigma^{2}}\).
Peters. The ecological implications of body size. Cambridge University Press, Cambridge, United Kingdom, 1983.
C. J. Pennycuick. Newton rules biology. Oxford University Press, Oxford United Kingdom, 1992.
A. Gelman and J. Hill. Data analysis using regression and multilievel / hierarchical modeling. Cambridge University Press, Cambridge, United Kingdom, 2009.
Exercise 3: Using for loops
Write a code fragment to set vague normal priors for 5 regression coefficients – dnorm(0, 10E-6) – stored in the vector b.
for(i in 1:5){
b[i] ~ dnorm(0, .000001)
}
Exercise 4: Coding the logistic regression
Write R code (algorithm 3) to run the JAGS model (algorithm 3) and estimate the parameters, \(r\), \(K\)\(\sigma\), and \(\tau\). We suggest you insert the JAGS model into this R script using the sink() command as shown in algorithm 4. because this model is small. You will find this a convenient way to keep all your code in the same R script. For larger models, you will be happier using a separate file for the JAGS code. Here is the joint distribution for our logistic model again, with the priors updated from exercise 2 and \(\tau\) expressed as a derived quantity:
We use the sink() command to create a JAGS script from our joint distribution. This file is created within R and saved in the working directory. Please note that the outer set of brackets are only required when running this code within an R markdown document (as we did to make this answer key). If you are running them in a plain R script, they are not needed.
{ # Extra bracket needed only for R markdown files
sink("LogisticJAGS.R")
cat("
model {
# priors
K ~ dunif(0, 4000)
r ~ dunif (0, 2)
sigma ~ dunif(0, 50)
tau <- 1/sigma^2
# likelihood
for(i in 1:n) {
mu[i] <- r - r/K * x[i]
y[i] ~ dnorm(mu[i], tau)
}
}
",fill = TRUE)
sink()
} # Extra bracket needed only for R markdown files
Then we run the remaining commands discussed in the JAGS Primer. Note that jm is calling the JAGS script LogisticJAGS.R we just created.
#Be sure to do this sorting. It will be important for plotting later.
Logistic <- Logistic[order(Logistic$PopulationSize),]
inits = list(
list(K = 1500, r = .2, sigma = .01),
list(K = 1000, r = .15, sigma = .5),
list(K = 900, r = .3, sigma = .01))
data = list(
n = nrow(Logistic),
x = as.double(Logistic$PopulationSize),
y = as.double(Logistic$GrowthRate))
n.adapt = 1000
n.update = 10000
n.iter = 10000
jm = jags.model("LogisticJAGS.R", data = data, inits = inits, n.chains = length(inits), n.adapt = n.adapt)
## mean sd 2.5% 50% 97.5%
## K 1.236753e+03 6.251288e+01 1128.5558515 1.231539e+03 1.374330e+03
## r 2.008537e-01 9.705856e-03 0.1816925 2.008424e-01 2.199169e-01
## sigma 2.864319e-02 3.049836e-03 0.0234551 2.838107e-02 3.522825e-02
## tau 1.259384e+03 2.608176e+02 805.7827817 1.241488e+03 1.817713e+03
## Rhat
## K 1
## r 1
## sigma 1
## tau 1
Exercise 5: Understanding coda ojects
Look at the first six rows of the data frame.
Find the maximum value of \(\sigma\).
Find the mean of r for the first 1000 iterations.
Find the mean of \(r\) after the first 1000 iterations.
Make two publication quality plots of the marginal posterior density of K, one as a smooth curve and the other as a histogram.
Compute the probability that K < 1600. Hint: what type of probability distribution would you use for this computation? Investigate the the dramatically useful R function ecdf().
Compute the probability that 1000 < K < 1300.
Compute the .025 and .975 quantiles of K. Hint–use the R quantile() function.
# Make a data frame
df = as.data.frame(rbind(zm[[1]], zm[[2]], zm[[3]]))
# Look at the first six rows:
head(df)
# Find the mean of r for the first 1000 iterations
mean(df$r[1: 1000])
## [1] 0.1999867
# Find the mean of r for the last 1000 iterations
nr = length(df$r)
mean(df$r[(nr - 1000): nr])
## [1] 0.2009966
# Make a publication quality plot of the marginal posterior distribution of K as a smooth curve.
# More iterations would produce a smoother cruve) and as a histogram.
plot(density(df$K), main = "", xlim = c(1000, 1500), xlab = "K")
# Find the probability that the parameter K exceeds 1600
1 - ecdf(df$K)(1600)
## [1] 6.666667e-05
# Find the probability that the parameter 1000 < K < 1300
ecdf(df$K)(1300) - ecdf(df$K)(1000)
## [1] 0.8500333
# Compute the .025 and .975 quantiles of K
quantile(df$K, c(.025, .975))
## 2.5% 97.5%
## 1128.556 1374.330
Exercise 6: Summarizing coda objects
Summarize the coda output from the logistic model with 4 significant digits. Include Rhat and effective sample size diagnostics (more about these soon). Summarize the coda output for \(r\) alone.
## mean sd 2.5% 50% 97.5% Rhat n.eff
## r 0.2009 0.0097 0.1817 0.2008 0.2199 1 7576
Exercise 7: Using MCMCchains
Extract the chains for r and sigma as a data frame using MCMCchains and compute their .025 and .975 quantiles from the extracted object. Display three significant digits. Make a publication quality histogram for the chain for sigma. Indicate the .025 and .975 Bayesian equal-tailed credible interval value with vertical red lines. Overlay the .95 highest posterior density interval with vertical lines in blue. This is a “learn on your own” problem intended to allow you to rehearse the most important goal of this class: being sufficiently confident to figure things out. Load the package HDinterval, enter HDInterval at the console and follow your nose. Also see Hobbs and Hooten Figure 8.2.1.
ex = as.data.frame(MCMCchains(zm, params = c("r", "sigma")))
class(ex)
Fig. 4. Median and 95% credible intervals for predicted growth rate and posterior density of K.
The intervals are exactly overlapping in this particular case. Such overlap will not always occur. Equal-tailed intervals based on quantiles will be broader than highest posterior density intervals when marginal posterior distributions are asymmetric. There is a 95% probability that the true value of the mean per-capita population growth rate falls between these lines. Not that this differs from the prediction of a new observation of population growth rate, which would have much broader credible intervals.
Exercise 9: Creating caterpillar plots with MCMCplot
Use the MCMCplot function to create caterpillar plots to visualize the posterior densities for \(r\), and \(\sigma\) using the coda object zm from earlier. Then use MCMCplot() to explore different plotting options. There are lots of these options, including ways to make the plots publication quality, show overlap with zero, etc.
MCMCplot(zm, params = "mu")
Exercise 10: Assessing convergence
Rerun the logistic model with n.adapt = 100. Then do the following:
Keep the next 500 iterations. Assess convergence visually with MCMCtrace() and with the Gelman-Rubin, Heidelberger and Welch, and Raftery diagnostics.
Update another 500 iterations and then keep 500 more iterations. Repeat your assessment of convergence.
Repeat steps 1 and 2 until you feel you have reached convergence.
Change the adapt phase to zero and repeat steps 1 – 4. What happens?
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1.01 1.01
## r 1.00 1.00
## sigma 1.00 1.01
## tau 1.00 1.01
##
## Multivariate psrf
##
## 1
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.952
## r passed 1 0.834
## sigma passed 1 0.255
## tau passed 1 0.347
##
## Halfwidth Mean Halfwidth
## test
## K passed 1246.319 1.52e+01
## r passed 0.199 1.97e-03
## sigma passed 0.029 3.92e-04
## tau passed 1233.763 3.26e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.467
## r passed 1 0.442
## sigma passed 1 0.465
## tau passed 1 0.374
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.36e+01
## r passed 2.00e-01 1.77e-03
## sigma passed 2.87e-02 3.27e-04
## tau passed 1.25e+03 2.72e+01
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.917
## r passed 1 0.772
## sigma passed 1 0.974
## tau passed 1 0.983
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.17e+01
## r passed 2.00e-01 1.26e-03
## sigma passed 2.88e-02 3.67e-04
## tau passed 1.24e+03 2.89e+01
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1 1
## r 1 1
## sigma 1 1
## tau 1 1
##
## Multivariate psrf
##
## 1
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.368
## r passed 1 0.431
## sigma passed 1 0.812
## tau passed 1 0.762
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 1.55e+01
## r passed 2.00e-01 2.07e-03
## sigma passed 2.86e-02 3.69e-04
## tau passed 1.26e+03 2.61e+01
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.920
## r passed 1 0.738
## sigma passed 1 0.666
## tau passed 1 0.664
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 12.60234
## r passed 2.00e-01 0.00180
## sigma passed 2.84e-02 0.00032
## tau passed 1.28e+03 28.41852
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.766
## r passed 1 0.817
## sigma passed 1 0.375
## tau passed 1 0.216
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.25e+03 2.30e+01
## r passed 2.00e-01 3.16e-03
## sigma passed 2.85e-02 3.07e-04
## tau passed 1.27e+03 2.57e+01
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## You need a sample size of at least 3746 with these values of q, r and s
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## K 1 1
## r 1 1
## sigma 1 1
## tau 1 1
##
## Multivariate psrf
##
## 1
heidel.diag(zm.short)
## [[1]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.636
## r passed 1 0.665
## sigma passed 1 0.941
## tau passed 1 0.684
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 4.622118
## r passed 2.01e-01 0.000616
## sigma passed 2.86e-02 0.000118
## tau passed 1.26e+03 9.142824
##
## [[2]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.189
## r passed 1 0.231
## sigma passed 1 0.411
## tau passed 1 0.485
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 4.865863
## r passed 2.01e-01 0.000618
## sigma passed 2.86e-02 0.000114
## tau passed 1.26e+03 9.497840
##
## [[3]]
##
## Stationarity start p-value
## test iteration
## K passed 1 0.0846
## r passed 1 0.0725
## sigma passed 1 0.3529
## tau passed 1 0.4007
##
## Halfwidth Mean Halfwidth
## test
## K passed 1.24e+03 5.393596
## r passed 2.01e-01 0.000633
## sigma passed 2.86e-02 0.000111
## tau passed 1.26e+03 8.870200
raftery.diag(zm.short)
## [[1]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 12 14032 3746 3.75
## r 12 11052 3746 2.95
## sigma 4 4792 3746 1.28
## tau 6 6878 3746 1.84
##
##
## [[2]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 6 7003 3746 1.87
## r 9 9497 3746 2.54
## sigma 4 4955 3746 1.32
## tau 5 5577 3746 1.49
##
##
## [[3]]
##
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95
##
## Burn-in Total Lower bound Dependence
## (M) (N) (Nmin) factor (I)
## K 10 10206 3746 2.72
## r 7 7820 3746 2.09
## sigma 4 4713 3746 1.26
## tau 5 6078 3746 1.62
ADVANCED Exercise 11: Coding the logistic regression to run in parallel
Append R code (algorithm 5) to the script you made in exercise 4 to run the JAGS model (algorithm 2) in parallel and estimate the parameters, \(r\), \(K\)\(\sigma\), and \(\tau\). Use the proc.time() function in R to compare the time required for the sequential and parallel JAGS run. If your computer has 3 threads, try running only 2 chains in parallel when doing this exercise. If you have fewer than 3 threads, work with a classmate that has at least 3 threads.
Looks as if the parallel model runs 2.63 times faster. This factor should increase the more iterations you run (why?).
You might be very tempted, in the name of reproducible science, to set the seed inside each worker to the same value. What happens when you do this? Compare the results from this run to the previous run using MCMCsummary() and MCMC(trace).
## mean sd 2.5% 50% 97.5%
## K 1.238468e+03 62.688701935 1133.0447273 1.232635e+03 1.381504e+03
## r 2.005799e-01 0.009655304 0.1815338 2.005211e-01 2.192294e-01
## sigma 2.865261e-02 0.003033771 0.0234590 2.835532e-02 3.536203e-02
## Rhat
## K 1
## r 1
## sigma 1
MCMCsummary(zmP)
## mean sd 2.5% 50% 97.5%
## K 1.239659e+03 6.227385e+01 1.134518e+03 1.233489e+03 1378.0986189
## r 2.005185e-01 9.584991e-03 1.817614e-01 2.004705e-01 0.2194044
## sigma 2.859709e-02 2.896285e-03 2.358206e-02 2.843383e-02 0.0348032
## tau 1.260124e+03 2.518494e+02 8.255847e+02 1.236885e+03 1798.1932583
## Rhat
## K NaN
## r NaN
## sigma NaN
## tau NaN
MCMCtrace(zm, pdf = FALSE)
MCMCtrace(zmP, pdf = FALSE)
Repeat this exercise with fixed initial values using algorithm 6.
